## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.4
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: lattice
## Loading required package: crayon
## 
## Attaching package: 'crayon'
## The following object is masked from 'package:ggplot2':
## 
##     %+%
## Loading required package: Biobase
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following object is masked from 'package:Matrix':
## 
##     which
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Parsed with column specification:
## cols(
##   Chr = col_character(),
##   Position = col_double(),
##   Marker = col_character(),
##   ref = col_character(),
##   alt = col_character()
## )

Manipulate Files

Read the Cross

##  --Read the following data:
##   571  individuals
##   30328  markers
##   10  phenotypes
## Warning in summary.cross(cross): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
##   (Perhaps it is in basepairs?)
##  --Cross type: f2
## Warning in summary.cross(crossobj): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
##   (Perhaps it is in basepairs?)
##     F2 intercross
## 
##     No. individuals:    571 
## 
##     No. phenotypes:     10 
##     Percent phenotyped: 100 94.7 94.7 87 99.8 99.8 99.8 96.5 96.5 96.5 
## 
##     No. chromosomes:    7 
##         Autosomes:      Fvb1 Fvb2 Fvb3 Fvb4 Fvb5 Fvb6 Fvb7 
## 
##     Total markers:      30328 
##     No. markers:        3890 4101 4745 4013 4451 5381 3747 
##     Percent genotyped:  100 
##     Genotypes (%):      AA:29.2  AB:35.6  BB:35.2  not BB:0.0  not AA:0.0

Optional: form linkage groups #{r, echo=F} newcrossobj <- formLinkageGroups(crossobj, max.rf=0.35, min.lod=6, reorgMarkers=TRUE, verbose=FALSE) newcrossobj <- orderMarkers(newcrossobj, use.ripple = TRUE, window = 5, map.function = c("kosambi"), verbose = T) plot(newcrossobj) summaryMap(newcrossobj) #

QTL Mapping and Linkage

## [1] 0.4960421
## [1] 0.3184846 0.9235690

##         ind1       ind2 prop_match
## 1 52QKSimTco R9KUKOjnBl      0.924
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 30 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.

##               chr   pos CWT_WEO CWT_SAL1 CWT_SAL2
## AX-184840405 Fvb1 12035    26.0     53.8     60.2
## AX-184687419 Fvb2 18860    25.8     46.1     52.7
## AX-184908616 Fvb3 22445    26.7     47.9     56.4
## AX-184169910 Fvb4  3765    27.8     46.7     54.5
## AX-184843747 Fvb5   295    27.4     55.2     65.8
## AX-184677344 Fvb6  2525    27.5     52.3     60.3
## AX-184202548 Fvb7  4820    28.1     57.4     65.6
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 30 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.

##               chr   pos CCT_WEO CCT_SAL1 CCT_SAL2
## AX-184346603 Fvb1  8420    13.0     14.4    16.48
## AX-184687419 Fvb2 18860    14.2     16.5    18.85
## AX-184908616 Fvb3 22445    11.0     10.1    13.93
## AX-184169910 Fvb4  3765    11.5     12.7    15.78
## AX-184170905 Fvb5  6970    11.4     12.6    11.26
## AX-89849903  Fvb6 10950    11.1      9.7     7.44
## AX-184164903 Fvb7 10865    11.9     13.1    12.10
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 74 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.
## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.

##               chr   pos CAWT_WEO CAWT_SAL1 CAWT_SAL2
## AX-184335455 Fvb1 14285     31.3      70.0      70.9
## AX-166511927 Fvb2  3805     28.6      84.9      92.9
## AX-184612678 Fvb3  9615     29.6      90.8      98.2
## AX-123367232 Fvb4 19085     29.4      50.6      48.3
## AX-184595061 Fvb5 20095     29.4      66.8      67.7
## AX-184615096 Fvb6 23015     30.2      86.9      92.3
## AX-184590330 Fvb7   180     27.1      81.7      96.9

GWAS

## [1] 1.02

WEO

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -188.263   6:47:24      1           0
##     2      -175.97   6:47:24      1           0
##     3      -173.329   6:47:25      2           0
##     4      -172.915   6:47:26      3           0
##     5      -172.892   6:47:27      4           0
##     6      -172.89   6:47:28      5           0
##     7      -172.89   6:47:28      5           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -222.416   6:50:28      1           0
##     2      -214.567   6:50:28      1           0
##     3      -211.774   6:50:29      2           0
##     4      -211.138   6:50:30      3           0
##     5      -211.076   6:50:30      3           0
##     6      -211.069   6:50:31      4           0
##     7      -211.069   6:50:31      4           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -164.067   6:53:21      1           0
##     2      -158.203   6:53:21      1           0
##     3      -157.561   6:53:22      2           0
##     4      -157.52   6:53:22      2           0
##     5      -157.519   6:53:23      3           0
## Performing GWAS evaluation

SAL1

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -138.928   6:55:49      0           0
##     2      -115.797   6:55:50      1           0
##     3      -114.631   6:55:51      2           0
##     4      -114.573   6:55:52      3           0
##     5      -114.573   6:55:53      4           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -213.835   6:58:32      0           0
##     2      -198.408   6:58:33      1           0
##     3      -193.384   6:58:34      2           0
##     4      -191.709   6:58:35      3           0
##     5      -191.323   6:58:35      3           0
##     6      -191.225   6:58:36      4           0
##     7      -191.2   6:58:37      5           0
##     8      -191.194   6:58:37      5           0
##     9      -191.192   6:58:38      6           0
##     10      -191.192   6:58:39      7           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -54.232   7:1:19      1           0
##     2      -26.0886   7:1:20      2           0
##     3      -24.3347   7:1:21      3           0
##     4      -24.087   7:1:21      3           0
##     5      -24.073   7:1:22      4           0
##     6      -24.0721   7:1:23      5           0
## Performing GWAS evaluation

SAL2

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -110.883   7:4:3      0           0
##     2      -86.6082   7:4:4      1           0
##     3      -85.8537   7:4:5      2           0
##     4      -85.8348   7:4:6      3           0
##     5      -85.8347   7:4:6      3           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -208.738   7:6:42      1           0
##     2      -193.975   7:6:42      1           0
##     3      -188.208   7:6:43      2           0
##     4      -186.261   7:6:43      2           0
##     5      -185.89   7:6:44      3           0
##     6      -185.822   7:6:45      4           0
##     7      -185.809   7:6:45      4           0
##     8      -185.807   7:6:46      5           0
##     9      -185.807   7:6:47      6           0
## Performing GWAS evaluation

## iteration    LogLik     wall    cpu(sec)   restrained
##     1      -25.7674   7:9:41      1           0
##     2      6.386   7:9:42      2           0
##     3      7.01255   7:9:42      2           0
##     4      7.06207   7:9:43      3           0
##     5      7.06264   7:9:44      4           0
## Performing GWAS evaluation

Checking the outcome

## [1] 1.001463

## [1] 1.037738

## [1] 1.095772

## [1] 1.04106

## [1] 0.9822432

## [1] 0.9875107

## [1] 1.177989

## [1] 1.096744

## [1] 1.086727

Mutivariate attempt #```{r, echo=FALSE} #CWT_multiple years miss10 <- pheno[is.na(pheno$CWT_SAL2),] MixCWT_yr <- sommer::GWAS(fixed = cbind(CWT_SAL1,CWT_SAL2)Year, random=vs(ID, Gu= KA), rcov=~units, data=pheno[which(pheno$ID %in% rownames(KA)),], n.PC = 0, min.MAF=0, M= GO[!(rownames(GO) %in% miss10$ID),], gTerm = “u:ID”)

us7 <- as.data.frame(t(UniCWT_SAL2\(scores)) us7\)Marker <- rownames(us7) MI2_CWT_SAL2 <- merge(markers,us7,by=“Marker”,all.x = TRUE) colnames(MI2_CWT_SAL2) <- c(“Marker”,“Chrom”,“Position”,“ref”,“alt”,“CWT_SAL2 beta”,“p.val”,“CWT_SAL2 Fstat”,“R2”,“R2s”) manhattan(MI2_CWT_SAL2, pch=20,cex=.5, PVCN = “color score”) ```